\n",
"\n",
"Note: \n",
"This lecture is based on Chapter 1 of [Fundamentals of Numerical Computation by Driscoll, Braun] which is avaliable online \n",
"
\n",
"\n",
"Note: \n",
"These notes contain some (but not all) of the content of the lecture. Please let me know if there are any additions that you feel are particularly important and let me know if you find any errors. \n",
"
"
]
},
{
"cell_type": "markdown",
"id": "e61992d9",
"metadata": {},
"source": [
"Relative error of the phrase \"2.2 pounds of jam is the same as a kilogram\":"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "3678bd7c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0020967851631676034"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 2.20462262;\n",
"y = 2.2;\n",
"\n",
"abs( x- y) / x"
]
},
{
"cell_type": "markdown",
"id": "59a0d3f3",
"metadata": {},
"source": [
"- Review of last week,\n",
"- [ChimeIn](https://chimein2.cla.umn.edu/chime/6801/folder/24220)\n",
"- Whiteboard I\n",
" - Standard/scientific form \n",
" - general floating point numbers \n",
" - underflow, overflow, machine epsilon\n",
" - Examples: half, single, double precision \n",
"- Half precison = Float16 (Julia)\n",
"- Whiteboard II \n",
" - round-off error, relative error in floating point representations and the distribution of floating-point numbers,\n",
" - Conditioning \n",
" - Stability "
]
},
{
"cell_type": "markdown",
"id": "c929e6d6",
"metadata": {},
"source": [
"## Floating point numbers \n",
"\n",
"(standard/scientific form, Example: Planck's constant)\n",
"\n",
"For fixed $\\ell, k$, we define\n",
"\n",
"\\begin{align}\n",
" \\mathbb F_{k\\ell} := \n",
" &\\Bigg\\{ \n",
" (-1)^s (1 + f) 2^n : f = (0.b_1b_2....b_k)_2, \\quad n = c - (2^{\\ell-1}-1), \\quad c \\in \\{1,\\dots,2^{\\ell}-2\\}\\nonumber\n",
" \\Bigg\\} \\nonumber\\\\\n",
" %\n",
" &\\cup \\{\\pm0, \\pm\\infty, \\text{NaN}\\} \\cup \\{ \\text{subnormal numbers} \\} \\\\ \n",
"\\end{align}\n",
"\n",
"Here, \n",
"* $k$ is the *binary precision*, \n",
"* $s\\in \\{0,1\\}$ is the sign bit, so $(-1)^s$ is the *sign*, \n",
"* $1 + f$ is the *significand*, \n",
"* $f$ is the *mantissa*, also called the *fraction*, \n",
"* $n$ is the *exponent*, belongs to $\\{ 2 - 2^{\\ell-1}, ...., 2^{\\ell-1}-1 \\}$, \n",
"* $c$ is the *characteristic*, \n",
"* $2^{\\ell-1}-1$ is the *bias*\n",
"\n",
"This is stored in a computer as the following $\\ell + k + 1$ bits:\n",
"\n",
"\\begin{align} \n",
" \\texttt{bitstring} \\quad = \\quad s \\quad | \\quad c \\text{ written in binary} \\quad | \\quad b_1b_2....b_k\n",
"\\end{align} \n",
"\n",
"***Special values.*** \n",
"\n",
"When $c = 0$, $\\texttt{bitstring}$ represents the number $(-1)^s f \\cdot 2^{-(2^{\\ell-1} - 1)}$. When $f = 0$, this encodes $\\pm 0$ and when $f \\not=0$ these are called *subnormal numbers*,\n",
"\n",
"When $c = 2^{\\ell}-1$, $\\texttt{bitstring}$ represents $\\pm \\infty$ when $f = 0$ and $\\text{NaN}$ when $f\\not=0$."
]
},
{
"cell_type": "markdown",
"id": "556ad41c",
"metadata": {},
"source": [
"
Example. \n",
"\n",
"Write down all the floating point numbers with $k = 2$ in the interval $[1,2)$. What about the interval $[2,4)$? $[2^n,2^{n+1})$? \n",
"\n",
"There are $2^k$ equally spaced floating point numbers in any interval $[2^n,2^{n+1})$ (as long as $\\ell$ is sufficiently large)\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "7ee9b76d",
"metadata": {},
"source": [
"*Maximum number* is $2^{2^{\\ell-1} -1}( 2 - 2^{-k} ) \\approx 2^{2^{\\ell-1}}$. If the result of a computation exceeds this value, it is stored as $+\\infty$ and we have an ***OVERFLOW*** error,\n",
"\n",
"*Minimum positive (normal) number* is $2^{2-2^{\\ell-1}}$. If the result of a computation is between this value and zero, it is stored as $0$ and we have an ***UNDERFLOW*** error,\n",
"\n",
"*Machine precision* is $\\epsilon_{\\text{mach}} = 2^{-k}$ "
]
},
{
"cell_type": "markdown",
"id": "f067fca8",
"metadata": {},
"source": [
"
Examples.
"
]
},
{
"cell_type": "markdown",
"id": "851a9fd4",
"metadata": {},
"source": [
"In the following, we will write $1|\\ell|k$ for the number of bits used to represent the sign ($1$), the exponent ($\\ell$), and the significand ($k$). \n",
"\n",
"In practice, we have: \n",
"``Float16`` = half precision = $1|5|10$, \n",
"``Float32`` = single precision = $1|8|23$, and \n",
"``Float64`` = double precision = $1|11|52$, "
]
},
{
"cell_type": "markdown",
"id": "bc268740",
"metadata": {},
"source": [
"## Floating point arithmetic\n",
"\n",
"We now consider double precision $\\mathbb F := \\mathbb F_{52,11}$ and define a function $\\mathrm{fl} : \\mathbb R \\to \\mathbb F$ that maps real numbers to (one of) the closest points in $\\mathbb F$. \n",
"\n",
"Suppose that $x \\in [2^n, 2^{n+1}]$ (and $\\ell$ is sufficiently large). The floats in the interval $[2^n, 2^{n+1})$ are of the form \n",
"$(1 + z \\cdot 2^{-k}) \\cdot 2^n$ for $z = 0,\\dots,2^{k}-1$. Therefore, consectutive floats are spaced apart by $2^{n-k}$ and since $|x| \\geq 2^n$, we have\n",
"\n",
"\\begin{align}\n",
" \\frac\n",
" {| x - \\mathrm{fl}(x) |}\n",
" {|x|} \n",
" \\leq \\frac12\n",
" \\frac\n",
" {2^{n-k}}\n",
" {|x|} \n",
" \\leq \\frac12 2^{-k} = \\frac12 \\epsilon_{\\rm mach}\n",
"\\end{align}\n",
"\n",
"**Subtractive cancellation** when adding or subtracting numbers so that the result is much smaller in magnitude. For example, we know $x = -1.0000000000000000000147$ to $23$ significant digits but $x + 1 = 1.47 \\times 10^{-22}$ to only $3$ significant digits.\n",
"\n",
"This concept can be illustrated using (relative) condition number: suppose we are interested in a quantity $f(x)$ for $x\\in \\mathbb R$. In practice, we can only store $\\widetilde{x} := \\mathrm{fl}(x)$ in a computer (recall that $\\widetilde{x} = x(1 + \\epsilon)$ for some $|\\epsilon|\\leq \\frac12 \\epsilon_{\\mathrm{mach}}$). Then the ratio of the relative error in the output to the relative error in the input in the limit as the precision goes to zero is known as the *(relative) condition number* $\\kappa_f(x)$: \n",
"\n",
"\\begin{align}\n",
"\\frac\n",
" {\\left|\\frac\n",
" {f(x) - f(\\widetilde{x})}\n",
" {f(x)} \\right|}\n",
" {\\left|\\frac\n",
" {x-\\widetilde{x}}\n",
" {x}\\right|}\n",
" %\n",
" &= \\left|\n",
" \\frac\n",
" {f(x) - f(x + \\epsilon x)}\n",
" {\\epsilon f(x)}\n",
" \\right| \\nonumber\\\\\n",
" %\n",
" &\\to \\kappa_f(x) \\quad \\text{as } \\epsilon \\to 0\n",
"\\end{align}\n",
"\n",
"That is, the (relative) condition number is approximately the constant of pr0portionality between the relative error in the output and the relative error of the input:\n",
"\n",
"\n",
"\\begin{align}\n",
" \\left|\\frac\n",
" {f(x) - f(\\widetilde{x})}\n",
" {f(x)} \\right|\n",
" \\approx \\kappa_f(x) |\\epsilon|\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"id": "a4a292e6",
"metadata": {},
"source": [
"
Example. \n",
"\n",
"The condition number of $f(x) = x + 1$ is $\\kappa_f(x) = \\left|\\frac{x}{x+1}\\right|$. This blows up when $x\\approx -1$ which explains subtractive cancellation that we mentioned above.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "55eab8b9",
"metadata": {},
"source": [
"If the condition number is large, we say the problem is *ill-conditioned* or the problem is badly conditioned\n"
]
},
{
"cell_type": "markdown",
"id": "5e6b33ff",
"metadata": {},
"source": [
"## ``Float16`` "
]
},
{
"cell_type": "markdown",
"id": "f502f414",
"metadata": {},
"source": [
"In Julia you can specify the precision by using ``FloatN()`` where $N = 16, 32,$ or $64$. e.g."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "862c6016",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.140625"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Float16( π );\n",
"x == 3.140625;\n",
"\n",
"(1 + 1/2 + 1/16 + 1/128)*2"
]
},
{
"cell_type": "markdown",
"id": "ea0ee94a",
"metadata": {},
"source": [
"You can see how this is represented in binary using ``bitstring()``:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "92002d0c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\"0100001001001000\""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bs = bitstring( x )"
]
},
{
"cell_type": "markdown",
"id": "1390956f",
"metadata": {},
"source": [
"Here, the first bit is the sign and the next 5 bits correspond to the characteristic:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "52d93a7c-adf9-4ec7-94f7-df326fa4a297",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(s, c) = (0, 16)\n"
]
}
],
"source": [
"s = parse(Int, bs[1:1], base=2);\n",
"c = parse(Int, bs[2:6], base=2);\n",
"@show ( s, c );"
]
},
{
"cell_type": "markdown",
"id": "430ed6d6",
"metadata": {},
"source": [
"The final bits ``bs[7:16]`` $:= b_1b_2 ... b_{10}$ represent the significand/mantissa in the following way \n",
"\n",
"\\begin{align}\n",
"f &= \\big(0 . b_1 b_2 b_3 b_4 b_5 b_6 b_7 b_8 b_9 b_{10}\\big)_2 \\nonumber\\\\\n",
"%\n",
"&= \\sum_{i=1}^{10} b_i 2^{-i}\n",
"\\end{align}\n",
"\n",
"That is, ``bs[7:16]`` is $2^{10} \\times$ larger than $f$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3935fccd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.5703125"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = 2^(-10) * parse(Int, bs[7:16], base=2)"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "6ef453a8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float16}:\n",
" 1.0\n",
" 1.57"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# exponent is c - 15 (characteristic minus bias)\n",
"# significand is 1 + f\n",
"[ exponent(x), significand(x) ]"
]
},
{
"cell_type": "code",
"execution_count": 68,
"id": "b755c179",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.140625"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y = (-1)^s * 2.0^( c - 15 ) * ( 1 + f )"
]
},
{
"cell_type": "markdown",
"id": "ac9fd6c0",
"metadata": {},
"source": [
"Recall that ``x = Float16( pi )``"
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "fbdac3e5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x == y"
]
},
{
"cell_type": "markdown",
"id": "a5999cd3",
"metadata": {},
"source": [
"Error in this representation: "
]
},
{
"cell_type": "code",
"execution_count": 71,
"id": "11245fb0",
"metadata": {},
"outputs": [],
"source": [
"using Plots \n",
"using LaTeXStrings "
]
},
{
"cell_type": "code",
"execution_count": 72,
"id": "c2f2dbc9",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"